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ABSTRACT 

Globular clusters should be born with significant numbers of stellar- mass black holes (BHs). It 
has been thought for two decades that very few of these BHs could be retained through the cluster 
lifetime. With masses ~ 10 Mq, BHs are ~ 20 times more massive than an average cluster star. They 
segregate into the cluster core, where they may eventually decouple from the remainder of the cluster. 
The small- N core then evaporates on a short timescale. This is the so-called Spitzer instability. Here 
we present the results of a full dynamical simulation of a globular cluster containing many stellar-mass 
BHs with a realistic mass spectrum. Our Monte Carlo simulation code includes detailed treatments of 
all relevant stellar evolution and dynamical processes. Our main finding is that old globular clusters 
could still contain many BHs at present. In our simulation, we find no evidence for the Spitzer 
instability. Instead, most of the BHs remain well-mixed with the rest of the cluster, with only the 
innermost few tens of BHs segregating significantly. Over the 12 Gyr evolution, fewer than half of the 
BHs are dynamically ejected through strong binary interactions in the cluster core. The presence of 
BHs leads to long-term heating of the cluster, ultimately producing a core radius on the high end of 
the distribution for Milky Way globular clusters (and those of other galaxies). A crude extrapolation 
from our model suggests that the BH-BH merger rate from globular clusters could be comparable to 
the rate in the field. 

Subject headings: binaries: close — globular clusters: general — Gravitational waves — Methods: 
numerical — Stars: kinematics and dynamics 



1. INTRODUCTION 

Typical globular clusters (GCs) should form ~ 100 — 
1000 BHs within ~ 3 Myr and could retain most of 
them initially, if t heir natal kicks are sufficiently low (see 
iWong et aljl2012l and references therein). With masses 
~ 10 M Q , these BHs become the most massive objects in 
the cluster after just ~ 10 Myr, so their dynamics will be 
very different than that of typical stars. The presence of 
BHs can affect the overall structur e and evolution of the 
parent cluster (|Mackev et al. 1 120081) . BHs accreting from 
a stellar companion can be visible as bright X-ray bi- 
naries (XRBs), which are in princi ple detectable both i n 
our own and other nearby galaxies (jKalogera eEaDHoM). 
Merging BH-BH binaries are key sources of gravitational 
waves (GW) which should be detectable by upcoming 
interf erometers such as Advanced LIGO (jHarrv et al.l 

It is well known that the formation rate per unit 
mass of XRBs is higher in cluster s by orders of mag - 
nitude compared to the field (e.g., iPoolev et al.l [2003). 
This indicates that dynamics must play an essen- 
tial role in producing cluster XRBs. Prior to 2007, 
however, there had not been a single detection of 
a black hole XRB inside a GC, although many had 
been identified in galactic fields. This appeared to 
agree with many theoretical studies suggesting that 
essentially all BHs within a star cluster should be 
ejected th rough dynamical interactions on a timescale 
- 10 9 yr (Kulkarni et al.l Il993t iSigurdsson k Hernquistl 



1993; 



2006; 



Porteeies Zwart k McMillanl 120001: lO'Learv et all 
Baneriee et al.l 120101 ). Key to all these previous 



m.morscher@u. northwcstern.edu 
s- umbreit @ northwestern . edu 
w-farr@northwestern.edu 
rasio@northwestern.edu 



studies is the expectation that BHs will segregate rapidly 
through dynamical friction, on a timescale ~ 100 Myr, 
and will succumb to the so-cal led Spitzer instability 
(| Spitz er 1969; Kulkarni et al.lll993l ). i.e., dynamically de- 
couple from the cluster by forming a central subclus- 
ter consisting primarily of BHs. The relaxation time for 
this small- N sub-cluster of BHs is very short, leading to 
prompt core collapse and evaporation. Through dynam- 
ical interactions, some BHs will be ejected in the form 
of tight BH-BH binaries that merge via GW emission 
within a Hubble time. 

This scenario was first proposed bas e d on sim- 
ple analytic estimates bv IKulkarni et ail (|1993f ) and 
ISigurdsson k Hernquistl ( 19931 ). The first direct N- 
body simulations of this effect used up to JV ~ 

10 4 particles (e.g iPortegies Zwart k McMillanl 120001: 
iMerritt et al.l 12004) . Larger direct A-body simulations 
bv lBaneriee et all (2010) used N ~ 10 5 , but with a single 
black hole mass (10 M s ) and no primordial b i naries . Us- 
ing simple dynamical models, lO'Learv et all (|2006| ) and 
iSadowski et al.l (|2008| ) studied the evolution of popula- 
tions of BHs that were either completely decoupled from 
or in equi librium with the clus ter, respectively (see dis- 
cussion m iDowning et~aTll2010D . 

Monte Carlo (MC) methods have made it possi- 
ble to model realistic GCs self- consistently with N ~ 

10 5 — 10 6 and sig nifica nt primordial binary fr actions 
(e.g.. iGiersz et all [2001 iChatteriee et alJ l2010h . The 
most realistic GC models with BHs to date are from 
IDowning efafl (|2010Ll2011l ). who used a Henon-type MC 
code to simulate clusters with N = 5 x 10 5 stars, a dis- 
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tribution of BH masses, and primordial binaries. These 
works used analytic cross sections to determine the re- 
sults of dynamical interactions, rather than direct inte- 
gration. 

Previous studies suggest that most BHs are ejected 
on a timescale ~ 10 yr. Hence, old GCs should 
have very few, if any, BHs left. However, in 2007 
the first candidate B H XRB inside a GC was de- 
tected in NGC 4472 (|Maccarone et all 120071) . Since 
then, several additional BH ca ndidates have been found 
in clusters in ot her galaxies ([Brassington et all 120101: 
Shih et all 120101: iBarnard et al.l 12011b iMaccarone et al.l 
20111) . iStrader et all 1)2012^ discovered two stellar-mass 



BHs in a Milky Way (MW) GC (M22). Assuming these 
BHs are accreting from white dwarf (WD) companions, 
and using cal c ulated formation and survi val rates from 
llvanova et all (2010), S trader et all (|2012l ) estimate that 
M22 has ~ 5 - 100 BHs. 

Furthermore, there have been a few recent theoreti- 
cal suggestions that significant nu mbers of BHs could 
still remain in some old clu sters dMackey et al. 120081 : 
IMoodv fc Sigurdssonl 12001 ) . IMackev et al.l (120081) used 
iV-body simulations with BHs to explain the radius-age 
trend in the clusters in the Magellanic Clouds; with dif- 
ferent initial retention fractions for BHs, they were able 
to reproduce the trend of increasing spread in core radius 
with age in these systems. In some models, they retained 
up to ~ 100 BHs over a Hubble time. 

Here we re-examine the BH retention question based 
on a realistic, large- N, fully self-consistent MC model. 
We find that that at least some old MW clusters may 
indeed retain a large fraction of their primordial BHs. 
This dramatically different picture for the fate of BHs in 
GCs may have implications for both BH XRBs and the 
production of merging BH-BH binaries. 

2. METHOD 

We use a Henon-type MC method to self-consistently 
model the evolution of star clusters due to the effects 
of two-body relaxation, strong binary scattering encoun- 
ters, stellar collisions, single and binary stellar evolu- 
tion, and mass loss from the Galactic tidal field. A 
detailed description of our code, as well as examples 
of its capabilitie s and comp a risons with o ther methods, 
can b e found inlJoshi et all (120001 12001D: iFregeau et all 
lf200l : lFregeau fc Rasiol (|2007D : IChatteriee et alJ ^OlO). 
The code has been well tested against direct iV-body 
models whenever possible. Since the dynamical evolu- 
tion of BHs in clusters is strongly dependent on inter- 
actions involving BH binaries, we perform direct calcu- 
lations of all strong 3-body (binary-single) and 4-body 
(binary-binar y) interactions using the small- N integra- 
tor Fewbody (jFregeau et al.l [2004D . These interactions 
are responsible for the hardening of BH-BH binaries and 
ejections of BHs from the cluster. Single star and bi- 
nary evolution are modeled using the routines of SSE 
and BSE (|Hurlev et al.ll2000l 120021 ) ■ Orbital energy loss 
from GW emission is handled within BSE for binaries 
retained in the cluster. For ejected systems, which are 
no longer evolved with our code, we use a si mplified 
times cale for GW inspiral in the weak field limit ()Petersl 
1 1964f ) . Neutron stars and BHs receive natal kicks with 
velocities drawn from a Maxwellian distribution with 
fi=265 km s . For BHs, the kick velocity is lowered 



according to the amount of material that falls back onto 
the final BH afte r the s upernova explosion, according to 
iBelczvnski et all (|2002D . 

We have recently added to our code a pre- 
scription for three-body binary formation, 
which is important for the evolution of BHs 
(Kul karni et al.l 119931: iSigurdsson fc Hernquistl 1993 : 
Portegies Zwart fc McMillanl 120001: lO'Learv et al.M2006 : 
Baneriee et al.l I2010D. W e follow a sim il ar pro cedure 



to llvanova et al. | ((20051) . llvanova et al.l (|2010D . and 
IQTearv et all (|2006f) to obtain an expression for the 
binary formation rate as a function of hardness ratio 



V = 



Gm\ 7712 

r p < m > a 2 



(1) 



We keep both the geometric and gravitational focus- 
ing contributions to the cross-section (in contrast to 
llvanova et alll2010l where the geometric part of the cross 
section for the third star to interact with stars 1 and 2 is 
dropped). For local number density, n, and average rel- 
ative velocity at infinity, Voo , the rate at which two stars 
(mi and 7772) form a binary with hardness 77 > 77 m i n 
through an interaction with a third star (7773) is given by 



r(?7 > 77„ 



= VSnVttf 

X (777! + 777 2 ) 5 77 m f n 5 (l + 27/ min ) 



1 + 2ry n 



777l + m-2 + Tl3 

mi + m 2 



(2) 



As we expect only dynamically hard binaries to survive 
(Hcggiel ["l975l ). we only consider the formation of hard 
binaries with 77 > 5 = 7y m i n . We allow three-body binary 
formation only for BHs. When forming a three-body bi- 
nary, we choose a value for 77 from a distribution accord- 
ing to the differential rate, dr/d7y, with lower limit 7y m ; n . 
The rest of the properties of the system are calculated 
from conservation of momentum and energy. 

We have checked that our MC prescription produces 
binaries at a rate that is in agreement with the analytic 
rate from Eq. ([2]). We ha ve also done a set of tests us- 
ing a direct N-body code (jFarr fc Bertschingerll2007t ) to 
check that our prescription produces hard binaries at the 
correct rate. Using one of our cluster snapshots, we in- 
tegrated our system of BHs for a short period of time 
with direct TV-body, and compared our binary forma- 
tion probability prediction to the actual binary forma- 
tion probability in the direct integration. We find good 
agreement with the direct iV-body trials, which gives us 
confidence that we are actually producing hard binaries 
at the correct rate. 

3. THE EVOLUTION OF A CLUSTER WITH 
STELLAR-MASS BLACK HOLES 

We present the results of a cluster model starting with 
N = 3 x 10 5 stars following a King profile with Wq = 5, 
half-mass radius = 2.44 pc, metallicity of Z = 0.001, 
and initial binary fr action fp = 0.1 . We choose our stel- 
lar masses from the iKroupal (|2001l ) initial mass function 
ranging from 0.1 - 100 M@. The cluster has initial total 
mass M to t = 2.03 x 10 5 M@ and half-mass relaxation time 
T r h ~ 6.5 x 10 8 yr. The central escape speed if] 31 km 

Our cluster is well below the veloci ty dispersion limit of ~ 40 
km s _1 found in Miller & Davics (2012), above which growth of a 
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Figure 1. Mass functions for BHs present in the cluster at 0.93 
Gyr (solid black line) and at 12 Gyr (dotted grey line). The pop- 
ulation at 0.93 Gyr is representative of the original population of 
initially retained BHs, as some (about 100) are ejected from the 
cluster at formation from natal kicks. Between 0.93—12 Gyr, the 
most massive (m bh > 20 Mq) are ejected preferentially over the 
lower mass BHs because they have higher interaction rates. 

s — 1 . Only about 12% of the BHs formed in the cluster re- 
ceived natal kick speeds above this value. We choose ou r 
remnant masses according to [Bclczvns ki et al.l ([2002), 
which produces BH masses in the range ^5 — 30 Mq 
for Z = 0.001. We form about 700 BHs in total, of 
which about 600 are retained initially (the remainder are 
ejected by natal kicks). The BH mass distribution at an 
early time is shown in Figure [T] 

Within a few Myr, the BHs begin to segregate, lead- 
ing to a central collapse by about 400 Myr. Formation 
of three-body binaries and their subsequent interactions 
lead to repeated core oscillations (see Figure [5]) . After 
about 300 Myr, strong binary interactions involving the 
mass-segregated BH population start to become dynam- 
ically important, and the rate of ejection of BHs (both 
single and binary) increases abruptly. Ejections continue 
through the end of the simulation, but the rate slows 
down over time. The evolution of the numbers of single 
and binary BHs retained in and ejected from the cluster 
is shown in Figure [3J For the entire simulation, most of 
the BHs are single; in fact, beyond about 300 Myr, there 
are typically no more than about 10 BH-binaries in the 
cluster. 

Statistics of the BHs at different times are shown in 
Table [T] By 12 Gyr, the cluster has ejected 202 single 
BHs, 33 BH-BH binaries, and 6 BH binaries with non- 
BH companions. Throughout the simulation, 13 BH-BH 
binaries merge due to GW emission; 6 of these merg- 
ers occur within the cluster, while the rest occur post- 
ejection. Most of the BH ejections and BH-BH mergers 
occur within about the first 6 Gyr of evolution. At 12 
Gyr, our model still has nearly 400 BHs, more than half 
of the initially-retained population. 

In Figure 2] we show the fractions of single BHs and all 
single stars in radial bins at several times. The BH frac- 
tion in the central bin, which always contains 20 BHs, 
grows to unity within about 600 Myr (left three panels), 
meaning that the innermost 20 objects are all BHs. Just 
outside the central bin, the BH fraction is typically less 
than 0.4, and it decreases to negligible fractions beyond 
about 1 pc. This indicates that, while the BHs do in- 
deed segregate to some extent, most of the BHs do not 
dynamically decouple from the cluster (i.e. they do not 
become Spitzer unstable). All but the innermost 20 or 

massive BH through successive mergers might occur. 



so most massive BHs remain well mixed with the cluster 
at all times. 

The most massive BHs tend to be preferentially ejected 
from the cluster (see Figure [Hand Table Q}. Nearly 75% 
of the ejected BHs have masses > 20 Mq, despite the fact 
that these more massive BHs are much less common than 
lower mass BHs. Since the most massive BHs sink the 
deepest, they tend to have the highest rates of strong 
interactions, which provide the energy needed to eject 
them from the cluster. 

We end at 12 Gyr, a typical age for MW GCs, with 
the cluster having N = 2.47 x 10 5 stars, M tot = 
1.05 x 1O 5 M , rh = 12.7 pc, and binary fraction f b = 
0.098. The final mass of our cluster is just slightly larger 
than the median value for MW GCs (M med w 8 x 10 4 ; 
iHeggie fc Hutll2003h . For our model we find an observa- 
tional core radius, r c ~ 5 — 7 pc, which falls within the 
high end of the core radius distribution of the MW GC 
system, and is also consistent with the range of core radii 
associated with old (~ 10 Gy r) GCs in the Mage llanic 
Clouds (see Figures 1 and 2 in IMackev et al.ll2008l) . 

4. DISCUSSION AND CONCLUSIONS 

The evolution and survivability of BHs in clusters, as 
well the effect that BHs have on their host cluster, will 
depend strongly on the degree to which the BHs are able 
to decouple from the cluster. Our MC method allows 
us to include realistic initial conditions as well as all the 
relevant physics for studying these types of systems in 
detail. 

In the most optimistic model of iMackev et al.l (|2008[ ) 
(run 4), about 50% of their BHs (s=s 100) are retained 
over ~ 10 Gyr. This is slightly less than our final reten- 
tion fr action (about 65 %) , but with N of three times 
that of Mackey et al. (2008), this amounts to more than 
a factor of three difference in the actual number o f BHs 
that we retain. In contrast with IMackev et al.l (2008) 
who found no BH-BH mergers within a Hubble time, we 
produce 13 mergers. In clusters with low central escape 
velocities, recoil kicks from strong dynamical encounters 
may tend to eject BH-BH binaries before they are tight 
enough to merge within a Hubble time. Although some 
of the Mackey et al. (2008) models do indeed have signifi- 
cantly lower escape velocities than the model we present, 
their run 4 actually has a comparable escape velocity, so 
this cannot reconcile the difference in merger rate. In- 
stead, the discrepancy may be explained by the larger 
number of BHs, as well as the inclusion of primordial bi- 
naries, resulting in a higher interaction rate in our sim- 
ulation, which is consistent with the larger number o f 
ejected BHs (but see discussion in iDowning etail pOlOl ) 
about the competing effects of hardening and destruction 
of BH-BH binaries, that go alon g with high BH i nterac - 
tion rates). We also compare to IDowning et a l. (2010), 
who track BH-BH mergers in a set of MC simulations. 
Their model 101ow75 is most similar to ours, with the 
same binary fraction and metallicity, and N = 5 x 10 5 , 
rh = 2 pc and T r h = 5.25 xlO 8 yr. They produced 6 
± 3 mergers (averaged over 10 simulations) within T#, 
about half as many as we produce in our simulation. 
Agreement to within a factor of two is reasonable, con- 
sidering their use of cross sections for predicting the out- 
comes of strong binary interactions (rather than direct 
integration) , which may overestimate the disruption rate 
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Figure 2. Evolution of the Lagrange radii of the BHs and all other objects. Plots show the radii enclosing 1%, 5%, 10%, and 50% of the 
mass, from bottom to top, for each category of object: BHs (solid curves) and non-BH (dotted curves). The left panel shows the evolution 
over the shorter interval corresponding to the grey band in the right panel. The BHs segregate from the rest of the cluster on a timescale 
of a few hundred Myr. The vertical grey dashed lines indicate the times when three-body binaries form. Interactions involving these 
hard binaries cause most of the oscillations of the innermost Lagrange radii. However, interactions involving hard binaries not created by 
three-body formation can also cause the same effect. 

clusters to date. This result is timely, considering the 
re cent discovery of two BH XRBs in a Milky Way GC 
bv lStrader et al.1 (|2012l ). who suggest that there may ac- 
tually be ~ 5 — 100 BHs in M22 at present. Our main 
conclusion is different from that of many other studies 
in the literature. This difference is not easily reconciled, 
but will be the subject of fut ure investigations. 

As has been suggested by iMackev et al.l (|2008t ) , the 
presence of BHs can indeed cause heating that can lead to 
significant core expansion, as we confirm with our model. 
The smaller cores observed in MW globulars may indi- 
cate l arger BH kicks than a ssumed in this work; intrigu- 
ingly, iRepetto et al.l (|2012f ) suggested such a change to 
the kick distribution on the basis of a population synthe- 
sis study of Galactic BHs. 
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Figure 3. Evolution of BH population retained in (top) and 
ejected from (bottom) the cluster. Each plot shows the number 
of single BHs (solid black line), BH-BH binaries (dashed red line), 
and BH- binaries with a non-BH companion (dotted blue line), ei- 
ther retained or ejected. The number of BH-BH binaries within 
the cluster tends to decrease over time until there are just a few. 
Both single BHs and BH-BH binaries are ejected efficiently from 
about 300 Myr until about 6 Gyr, at which point the ejection rate 
slows down significantly. This occurs in conjunction with a drop in 
the overall binary interaction rate, which is caused by the ongoing 
cluster expansion due to heating by the BHs. BH-other binaries 
tend to increase in number very slowly over time. At 12 Gyr, the 
cluster still has 395 BHs, more than half of the cluster's initial 
population. The majority of the BHs are single at all times. 
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for tight BH-BH binaries (jDowning et al.ll2010D . 

A crude extrapolation from our model can be used to 
estimate the rate of BH-BH mergers in a Milky Way 
equivalent galaxy (MWEG). In our model, the total 
merger rate is ~ 1 per Gyr. Our Galaxy may have had 
~ 300 GCs (about half of which have since dissolved). 
We therefore estimate a merger rate of ~ 0.3 per MWEG 
per Myr from star clusters. This is exceeds the estimated 
merger rate from pr imordial binaries in the galactic field 
(jAbadie et al.l r2010). Thus, our model indicates that it is 
important to include GCs in calculations of the BH-BH 
merger rate of the Universe. 

Our results indicate that at least some old GCs could 
have hundreds of stellar- mass BHs at present. Since 
nearly all of our BHs are single, our prediction is con- 
sistent with the small number of BH XRBs detected in 
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Table 1 

Properties of BH population at different evolutionary stages: 0.93 Gyr, 3.25 Gyr, 6.5 Gyr, 9.77 Gyr and 
12 Gyr. The table shows the number of single BHs (AT b bh)i number of BH-BH binaries (A^bh— bh)j 
number of BH-other (non-compact) binaries (TVbh— other); number of BHs with masses above 20 Mq 
(JVjH(m > 20 Mq)), average individual BH mass (m ave bh); that are retained in/ejected from the 
cluster, at the different times. We also show the number of BH-BH mergers (-/V mgr ) that have occurred 
up to the time given, either inside the cluster or post ejection. 
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Figure 4. The fraction of single BHs (solid black line) and all single stars (dashed red line) at several times; the remainder of the 
objects are binaries. The innermost bin contains 20 BHs, and the number of BHs inside each subsequent bin doubles (40, 80, etc.). The 
BH fraction in the central bin (containing 20 BHs) reaches unity by about 600 Myr (left panels), and then fluctuates between 0.4-1 for the 
rest of the simulation (right panels). Beyond the first bin, the BH fraction decreases, reaching negligible fractions beyond about 1 pc. The 
vertical grey dashed line shows the extent of the innermost 140 BHs (20 + 40 + 80 = 140 contained within the first three bins), which is at 
all times less than half of the retained BH population. 
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